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RAPID DYNAMICAL CHAOS IN AN EXOPLANETARY SYSTEM 
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ABSTRACT 

We report on the long-term dynamical evolution of the two-planet Kepler-36 system, which we studied 
through numerical integrations of initial conditions that are consistent with observations of the system. The 
orbits are chaotic with a Lyapunov time of only ~10 years. The chaos is a consequence of a particular set of 
orbital resonances, with the inner planet orbiting 34 times for every 29 orbits of the outer planet. The rapidity of 
the chaos is due to the interaction of the 29:34 resonance with the nearby first order 6:7 resonance, in contrast to 
the usual case in which secular terms in the Hamiltonian play a dominant role. Only one contiguous region of 
phase space, accounting for ~ 4.5% of the sample of initial conditions studied, corresponds to planetary orbits 
that do not show large scale orbital instabilities on the timescale of our integrations (~ 200 million years). The 
long-lived subset of the allowed initial conditions are those that satisfy the Hill stability criterion by the largest 
margin. Any successful theory for the formation of this system will need to account for why its current state is 
so close to unstable regions of phase space. 

Subject headings: celestial mechanics — planets and satellites: dynamical evolution and stability 



1. INTRODUCTION 

Despite the seeming regularity of the Solar system, the 
planetary orbits are kno wn to be chaotic with a Lyapunov 
time o f ^5 million years (lLaskar 1119891 [Wisdom & Sussmanl 
119921) . The hallmark of a chaotic system is sensitive depen- 
dence on initial conditions: two trajectories that start arbitrar- 
ily close to each other will diverge exponentially on a time 
scale known as the Lyapunov time. Chaos is seen not only in 
the orbits of the planets, but also among various satellites and 
minor bodies of the Solar Syste m (Wisdom 1985; Leca r et aTl 
1200 U iGoldreich & Rappaportl 120021) . However, among the 
hundreds of multiplanet systems known to exist around other 
stars, few have orbits measured precisely enough to determine 
definitively if chaos is present. There is evi dence that both the 
Kepler- 1 1 and 55Cnc sys tems are chaotic dGayon et al.ll2008l 
iMigaszewsk i et al. 2012), but in the case of 55Cnc, where the 
masses and inclinations are not well constrained, this conclu- 
sion is less certain. 

The Kepler— 36 system consists of a subgiant star of solar 
mass and the two transiting planets Kepler-36b and c, with or- 
bital periods of 13.8 and 16.2 d and mas ses of 4.1 and 7.5M®, 
respectivelv dCarter & Agol et al.ll2012b . All of the necessary 
parameters for integration — the bodies' positions and veloci- 
ties at a reference epoch, as well as their masses — have been 
measured precisely. This allows study of the true dynamical 
evolution of the system at a level of detail that has only been 
achieved for the Solar System and a handful of exopla netary 
systems dGozdziew ski 2005; Migasze wski et alj|2012l) . 

This Letter is organized as follows. Section [2]describes our 
integration methods and the set of initial conditions we use. 
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Our results on the chaotic behavior of the planets, including 
an explanation of its salient features, are given in Section [3] 
A stability analysis of the system is presented in Section @] 
We discuss briefly in Section |5]how our results may constrain 
models of formation of the system. 

2. NUMERICAL METHODS 

The effect on the transit light curves of dynamical interac- 
tions between the planets was modeled in a Bayesian fashion 
using prior knowledge of the host star obtained through astro- 
seismic analysis. This resulted in a posterior joint probabil- 
ity distribution for the bodies' masses an d the planetary posi- 
tions and velocities at a reference epoch (iCarter & Agol et all 
120121) . Initial conditions and masses drawn from this dis- 
tribution form a representative set because they sample the 
possible configurations of the planets consistent with the data 
in a statistically appropriate way. All of the correlations be- 
tween the uncertainties in these orbital parameters are natu- 
rally taken into account. 

We studied the dynamical evolution of the Kepler-36 
system through numerical integration of 10 4 initial condi- 
tions and massesQ Our prima ry integration scheme was a 
symplectic n-body mapping dWisdom & Holmanl 1199 ll) . 
We implemented Chambers' symple ctic corrector to im 
prove the accuracy of t h e inte gr ator dWisdom et al. 1 1 1996 
Chambers & Migliorini I 119971: lLaskar & Robutel I |2001 



Wisdom I I2006I) . resulting in a relative energy error over the 



course of our integrations of < 10~ 12 . 

"Stepsize ch aos" can be a source of error in symplec- 
tic integrators dWisdom & Holman I [19921 but it is negligi- 
ble when the timestep At is smaller than the shortest phys- 
ical timescale in the syst em by at least a factor of 10 or 20 
dRauch & Holmanl lT999). Each integration was carried out 
using a fixed timestep of At < l%Pb- 

Only Newtonian gravitational forces were considered; gen- 
eral relativity is unimportant due to the long timescale for 
relativistic precession (~ 10 8 days) relative to the secular 

6 These initial conditions are published with Carter & Agol et al. (2012). 
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FIG. 1. — Distribution of estimated Lyapunov times. Top: The distribu- 
tion for the long-lived initial conditions after 5 X 10 7 days (green), 10 7 days 
(black), 10 6 days (blue), and 10 5 days (red). Middle: The distribution for 
the entire set of initial conditions, same color scheme. Bottom: The trajec- 
tories with the fastest chaos correspond to orbits with higher initial mutual 
inclinations. 
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10 4 days. Tidal effects are discussed in Sec- 



3. RAPID DYNAMICAL CHAOS 

A direct way of determining if a particular trajectory is 
chaotic is to measure how perturbations of that trajectory 
grow in time. In a chaotic system, a perturbation will grow 
exponentially in the limit of t — > oo, with an e-folding 
time equal to the Lyapunov time. There are as many dis- 
tinct Lyapunov times as degrees of freedom of the system 
dLichtenberg & Lieberman 1 1 1992b . but in practice one usu- 
ally measures the shortest. In addition to a main trajectory, 
we integrate the tangent equations of the symplectic mapping 
which govern the behavior of infinitesimally small perturba- 
tions of that trajectory. 

The tangent equations were integrated forward by 10 7 days 
for each trajectory. Over 99% of these initial conditions, in- 
cluding the orbit that best fits the data, are chaotic. The rest do 
not show exponential divergence on the timescale of the inte- 
grations and therefore may be quasiperiodic. Figure Q~]shows 
the distribution of Lyapunov times for the entire set of initial 
conditions and for those belonging to the long-lived region 
(where we predict the true orbits must lie, see Section^). The 
typical Lyapunov time is several thousand days. A more per- 



tinent measure is the ratio of the Lyapunov time to the short- 
est orbital period. This ratio is ^300 for Kepler-36. In con- 
trast, the Lyapunov time of the Solar System is approximately 
2 x 10 7 orbits of Mercury. One of the most rapid Lyapunov 
times observed within the Solar System, that of the Saturnian 
moons Prometheus and Pan dora, is 2,000 Promethean orbits 
dFarmer & Goldreich ||2006|) . 

We confirmed that similar estimates of the Lyapunov time 
were obtained when using a smaller timest ep (At = 0.01 d) 
and w hen using a Bulirsch-Stoer integrator dBulirsch & Stoer I 
119911) . As expected, the distribution of Lyapunov times did 
not change when general relativity was included, which we 
mimicked using a dipole potential which p roduces the correct 
precession rate (iNobili & Roxburgh |[T98l . 

Figure Q] shows that the estimated Lyapunov times range 
over two orders of magnitude, with two clusters centered at 
approximately 300 and 4,000 days. This was unexpected, 
since chaotic zones in phase space are typically chara cterized 
by a single value of the minimum Lyapunov time dHenon I 
1983). Moreover, the estimated Lyapunov times change as the 
total integration time is increased, with a growing population 
centered on 300 days. This suggests that the initial conditions 
span two nearly disconnected regions of the chaotic zone in 
phase space, characterized by different estimates of the local 
Lyapunov time (which is not the absolute shortest Lyapunov 
time, in the infinite limit), and that trajectories are moving be- 
tween them. This will be returned to in Section |4] we turn 
now to understanding the rapidity of the chaotic behavior. 

3.1. Origin of the chaos 

To gain an understanding of the chaos associated with a 
Lyapunov time of several thousand days, we sought evidence 
for resonant behavior. Mean motion resonances (MMRs) be- 
come important when the ratio of the planetary mean motions 
is close to a rational number, rib/n c ~ (j + k)/j, where j 
and k are mutually prime integers. Numerically we find that 
the angles characterizing the 29:34 eccentricity-type MMR, 
9i = 34A c -29A h -5n7t and6> 2 = 34A c -29A h -5n7 c (and the 
linear combinations of 6>i and 82 appearing in the interaction 
Hamiltonian), to be behaving chaotically. Here A and m refer 
to the mean longitude and longitude of periastron of the or- 
bits. Chaotic alternations of these angles between circulation 
and libration are shown for a particular trajectory in Figure[2] 
It may be surprising that such a high order (large fc) MMR 
would be important, as resonance widths in phase space scale 
as k factors of eccentricities and inclinations (and therefore 
are usually negligible for nearly circular and coplanar orbits). 
However, the Laplace coefficients, which also factor into the 
widths, are large when the semimajor axis ratio of the plan- 
ets is near unity, as is the case for Kepler-36. This allows the 
29:34 resonance to be important. 

To understand why the 29:34 resonant angles should be- 
have chaotically, we appeal to the resonance overlap crite- 
rion dChirikov 1 119791) . This states that chaos will ensue if 
the angles 82 and 8\, when neglecting the interaction between 
them, are both analytically calculated to be librating in the 
same region of phase space. This can occur if the resonance 
widths, which are functions of the eccentricities, become large 
enough or if the resonance centers, determined to be where 
6\ = and 82 = 0, coincide. We find that the vast majority 
of the initial conditions exhibit oscillations with a period of 
several thousand days in the angle Wb — w c , indicating that 
w\y ~ w c , i.e. that the resonant islands are overlapping. 
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FIG. 2. — Chaotic evolution of the resonant angles 6\ = 34A C — 29A;, — 
5t*7(, and 82 = 34A C — 29A;, — 5ro c for a randomly chosen trajectory from 
the long lived region. The red overlaid points show a smoothed version of the 
black points to guide the eye. 

The approximate equations of motion for resonant angles 
resemble those of coupled pendulums driven by the oscilla- 
tions in the eccentricities and the periastra of the two planets. 
This correspondence can be used to show that the resulting 
chaotic behavior should hav e a Lyapunov time similar to the 
period of these oscillations (iHolman & Murray] 1 1996I) . For 
Kepler-36 this period is several thousand days, similar to the 
Lyapunov time, supporting this explanation of the origin of 
the chaos. 

This explanation alone is not surprising, as the connection 
between chaotic behavior and MMR overlap is well known. 
Why then is the Lyapunov time (relative to the smallest or- 
bital period) so short for this sytem? In previously known 
examples MMR overlap, the periodic driving of the resonant 
angles by the eccentricity and periastra are caused by secu- 
lar effects. In the Kepler-36 system, the driving is dominated 
by the effects of a nearby first-order MMR (Pb/P c ~ 6/7), 
which appear at a lower order (in the eccentricities) in the in- 
teraction potential. It can be shown that the proximity to the 
6:7 resonance does not alter the qualitative character of the 
standard secular solution for the eccentricities and periastra 
but that it does affect the frequencies of the oscillations sig- 
nificantly. Following Malhotra et al. (1989), we estimate the 
shorter of the two modified secular timescales to be ~ 8, 000 
days. We confirmed that the 6:7 terms are crucial for pro- 
ducing the correct behavior of the eccentricities and periastra 
numerically as well. 

The 6:7 inclination-type terms appear in the Hamiltonian 
at the same order in the inclinations as the secular terms and 
hence should not be as crucial as the 6:7 eccentricity terms. 
Trajectories that have been projected onto the invariable plane 
remain chaotic with a Lyapunov time of several thousand 
days, confirming that a coplanar model is sufficient to explain 
the chaos. 

A much smaller subset of the initial conditions shows 
chaotic behavior of the angles characterizing the 23:27 
eccentricity-type MMR. Figure[3]delineates the chaotic zones 
(Lyapunov time < 10 6 d) associated with separate MMR, for 
different eccentricities, demonstrating that both the 23:27 and 
29:34 resonances can be important. Moreover, the chaotic re- 
gion is not always confined to a small range in period ratio, 
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FIG. 3 . — Structure of the phase space near the Kepler-36 initial conditions. 
Red indicates the span of the long-lived trajectories on this plane and blue 
the entire sample. Black points correspond to trajectories that do not fit the 
data. Green lines indicate the nominal location of the 23:27 and 34:39 MMR. 
Only a b varies within each panel. Top: (e b ,e c ) ~ (0.02,0.0). Middle: 
(e&,e c ) ~ (0.021,0.006). Bottom: (e b ,e c ) ~ (0.036,0.0). 

indicating that the orbits of the planets may not always be sta- 
ble, or constrained. 

4. STABILITY OF THE KEPLER-36 SYSTEM 

There are two primary modes of instability in planetary 
systems. First, planets may have close encounters or colli- 
sions. An analytic criterion is available for two-planet sys- 
tems, based on conservation of angular mo mentum L and en- 
ergy E , which can exclude this possibility dMarchal & Bozisl 
119821) . in which case the planets are labeled Hill stable. The 
(sufficient but not necessary) criterion for stability is 

h = -cL 2 E > h C rit (1) 

where h cr u = 1 + 3 4 / 3 m(,m c TO„ 2 ^ 3 /(m b + ?ti c ) 4 / 3 + ... and 
c = 2(m* + to/, + nib) / (G 2 (mi,m c + Tn c m* + mjm,) 3 ). 
Here refers to the mass of the star. 

Even if the planets never have close encounters, repeated 
weak interactions can lead to a second type of instability, in 
which the gradual exchange of angular momentum and energy 
between the planets results in drastic orbital variations. This 
is known as Lagrange instability. There is no known analytic 
criterion for Lagrange stability, but numerical integrations can 
demonstrate that a given trajectory does not show unstable 
behavior on the timescale of the integration and therefore can 
be considered to be "long-lived". 

A preliminary stability analysis of the Kepler-36 system de- 
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termined that ~ 9% of th e initial conditions did no t satisfy 
the Hill stability criterion (ICarter & Agol et al.ll2012b . It was 
left to future work to determine whether the initial conditions 
were Lagrange long-lived. 

4.1. Lagrange stability analysis 

We integrated all 10 4 trajectories for 2.5 x 10 9 d and found 
that while only 30% showed Lagrange unstable behavior dur- 
ing the first ~ 10 8 d, this percentage had increased to ~ 75% 
by the end of the integrations. We classified Lagrange in- 
stability as variations in semimajor axis, eccentricity, and in- 
clination that were different from their running average by 
greater than 10%. In general, the trajectories that were La- 
grange unstable and those that did not show instability dur- 
ing these integrations were well mixed, but we could iden- 
tify a single contiguous region in phase space, accounting for 
~ 4.5% of the initial conditions, which had no unstable tra- 
jectories on these timescales. This region is characterized by 
lower eccentricities and inclinations, shows no systematic dif- 
ference in the goodness-of-fit to the data, and contains five 
of the candidate quasiperiodic trajectories. One hundred ran- 
domly chosen initial conditions from this "long-lived core" 
were integrated for 200 million years ( > 5 x 10 9 Pb). Of these, 
only four exhibited instability; all came from the borders of 
the long-lived core. 

The initial conditions belonging to the long-lived core sat- 
isfy 

h > h crlt + e (2) 

where the value of e ~ 0.0007 is determined numerically. 
In other words, orbits that are Lagrange long-lived must sat- 
isfy a criterion nearly identical to the Hill criterion, but with 
a slightly different critical value . This same relationship wa s 
found for Jupiter-mass planets dBarnes & Greenberg 1120061) . 
Our work suggests that the same link between Hill and La- 
grange stability applies across orders of magnitude in the pa- 
rameter m p i ane t I m st ar, though the value of e depends on this 
parameter. Contours of constant <Il = h— h cr u are indicated 
in Figure |4] 

The initial conditions belonging to the long-lived region 
are those with the smallest values of the angular momentum 
deficit (AMD), a quantity which parametrizes how far the or- 
bits are from purely circular and coplanar. It is straightfor- 
ward to show that is a function of the AMD for nearly 
circular, coplanar orbits, explaining the observed dependence. 

To explore the connection between and Lagrange insta- 
bility, we focused on a group of 475 Hill stable trajectories 
with a uniform distribution in dr,. There is a sharp transition 
at dL = 0.0007, where the trajectories change from being 
mostly Lagrange unstable to being entirely Lagrange long- 
lived. A smaller value of dL correlates with a shorter time 
to show Lagrange instability and with smaller minimum ap- 
proach distances between the two planets. 

It is not obvious what longer integrations would reveal re- 
garding the trajectories that appear to be long-lived on the 
timescale of our integrations but do not satisfy equation (O 
(i.e. are not in the long-lived core). Previous work indicates 
that unless these these orbits are protected by some resonance 
mechanism they too w ill show Lagrange instability on rela- 
tively short timescales dGladman Hl993tlBarnes & Greenberg I 
120071) . We predict that future observations will show that the 
true orbits of the Kepler— 36 planets belong to the long-lived 
core. 

By restricting the initial conditions to the long-lived 
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FIG. 4. — Relationship between Hill and Lagrange stability. Upper Plots: 
The region of orbital element space where orbits can be long-lived. Lagrange 
unstable orbits are shown in red, orbits that have not undergone Lagrange 
stability in 2.5 X 10 9 days are in black. Lower Plots: By grouping initial 
conditions based on their values of dr, we identify the Lagrange long-lived 



core we refine the system parameters. This results in 
median values and 84.2% and 15.8% confidence inter- 
vals for the masses and radii of the planets of m& = 



4.32Af ff 



7.84M 



1.49iL 



+0.035 



+019 „ _ 7 QAM + 033 R — 
3-0.20' m c — ' - 84JW ©-0.36' U b — -0.035' 

and R c = 3.68i?0lg ols- The ^ est ^ vames are m b = 
4.11M©,mc = 7.46Me,i? b = 1A6R®, and R c = 3.59i? e . 
The mutual inclination is constrained to be less than one de- 
gree. 

Although these planets seem to be on the brink of instabil- 
ity, the orbits are stable to small perturbations. In particular, 
we confirmed that there are hypothetical third planets (with a 
period of 2.6P C and M < 5M ffi ) that do not disrupt the two 
planet system on seven million year timescales. 

4.2. Transition to Lagrange instability 

We now return to the evolution of the distribution of Lya- 
punov times shown in Figure Q] We find that the onset 
of Lagrange unstable behavior occurs when the trajectory 
moves between the two peaks of the distribution (at ~ 300 
and ~4,000 days). A typical Lagrange unstable trajectory is 
shown in Figure [5] The eccentricities and semimajor axes 
remain nearly constant for about 10 5 days, after which the 
orbital elements vary dramatically. At the same time, the es- 
timated Lyapunov time changes from roughly 15,000 days to 
180 days. 
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FIG. 5. — Period ratio and eccentricity evolution (e;, in black, e c in red) 
for a Lagrange unstable orbit. When the orbit begins to change qualitatively 
in an erratic way, the slope of log(ci) changes, indicating a change in the 
estimated Lyapunov time, log(d) is the natural logarithm of the norm of 
distance between two initially nearby trajectories in phase space. 



The ability for trajectories to cross between the two nearly 
disconnected regions of the chaotic zone (which are char- 
acterized by different estimates of the local Lyapunov time 
and also are chaotic for different reasons ) can be under- 
stood as a consequence of r esonance overlap dMardling 120081 
iMurray & Holmaril 119971) . As Figure [3] shows, separate 
MMRs merge at higher eccentricities. Our hypothesis is that 
chaotic diffusion of eccentricities to higher values results in 
MMR overlap, at which point a bottleneck between the two 
regions is formed. Once this occurs, trajectories can enter the 



Lagrange unstable region where they explore a broader range 
of period ratios. This chaotic diffusion may occur for all of 
the chaotic trajectories, in which case the distribution of Lya- 
punov times will approach an asymptotic form with all trajec- 
tories belonging to the 300 day peak as longer integrations are 
performed. 

5. DISCUSSION 

We have presented a dynamical analysis of the Kepler-36 
system which has yielded several surprising results. The or- 
bits are chaotic with an extremely short Lyapunov time, yet 
some still manage to be long-lived. The closeness of this sys- 
tem to instability is an intriguing feature of Kepler-36. In par- 
ticular, did the planets form with orbits contained in the long- 
lived core? Alternatively, did some dissipative process drive 
them into this long-lived configuration? 

It would seem that tidal dissipation is negligible for this 
system, despite the proximity of the planets to the star. If tidal 
dissipation were important then the planets would have been 
on more eccentric orbits in the past, and, assuming the orbital 
separation was not very different at that time, it is unlikely 
such a configuration would have been stable. 

Alternatively, the planets could have been closer together 
if they were protected by a resonance. A commonly dis- 
cussed mechanism for forming compact multiplanet systems 
is convergent migration in the gase ous protoplanetary disk 
or the disk of remnant planetesimals (iTerquem & Papaloizou I 
12007b . Perhaps the Kepler-36 planets formed at large orbital 
distances, and migrated until being locked into a 6:7 res- 
onance. After migration ended, subsequent tidal evolution 
could have driven the planets apart to their current configu- 
ration. However, we recognize that it may be challenging to 
create closely packed r esonant systems th rough conventional 
migration mechanisms dRein et al. 1120 121) . though other stud- 
ies find that it is possible to produce systems like Kepler-3 6 



eph 

(llda & Linll2TM^Ogmara"et^l20Tof IPierens et al. N20TT1) . 
Therefore this system may contain important clues about the 
relevance of convergent migration to the formation and dy- 
namical evolution of planetary systems. 

We thank J. Wisdom for sharing with us his code that cal- 
culates coefficients in the disturbing function, as well as other 
members of the Kepler TTV team for helpful comments on 
the text. EA acknowledges support for this work was pro- 
vided by NSF Career grant AST-0645416. 
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